#---- Supplementary replication script for  ----

#Turnbull-Dugarte, López-Ortega & Hunklinger (2024)
#"Do citizens stereotype Muslims as the Illiberal 'Bogeyman'? Evidence from a Double-List Experiment"

# Packages
remove(list=ls())

#install.packages("pacman")
library(pacman)

p_load(here, tidyverse, haven, modelsummary, jtools,
       list, ggpubr, survey, xtable, interactions, viridisLite,
       ggdist, patchwork, estimatr, margins, grf, starbility, miceadds)



colors <- viridis(option = "plasma", 
                  begin = 0.2, 
                  end = 0.8, 
                  direction = -1, 
                  n = 4)


##Loading data##
list_dataNL <- read_dta("data/listNL.dta")
list_dataNL <-as.data.frame(list_dataNL)
list_dataNL<- list_dataNL  %>% 
  mutate(treatfac=as.factor(treated),
         IDvar=as.factor(IDvar),
         homonat=as.factor(homonat),
         round.fac=as.factor(round))
listNL_A <- list_dataNL %>% 
  filter(round==1)
listNL_B <- list_dataNL %>% 
  filter(round==2)
list_dataNL <- subset(list_dataNL, !(is.na(weight)))

list_dataDE <- read_dta("data/listDE.dta")
list_dataDE <-as.data.frame(list_dataDE)
list_dataDE<- list_dataDE  %>% 
  mutate(treatfac=as.factor(treated),
         IDvar=as.factor(IDvar),
         homonat=as.factor(homonat),
         round.fac=as.factor(round))
listDE_A <- list_dataDE %>% 
  filter(round==1)
listDE_B <- list_dataDE %>% 
  filter(round==2)
list_dataDE <- subset(list_dataDE, !(is.na(weight)))

list_dataUK <- read_dta("data/listUK.dta")
list_dataUK <-as.data.frame(list_dataUK)
list_dataUK<- list_dataUK  %>% 
  mutate(treatfac=as.factor(treated),
         IDvar=as.factor(IDvar),
         homonat=as.factor(homonat),
         round.fac=as.factor(round))
listUK_A <- list_dataUK %>% 
  filter(round==1)
listUK_B <- list_dataUK %>% 
  filter(round==2)

list_dataUS <- read_dta("data/listUS.dta")
list_dataUS <-as.data.frame(list_dataUS)
list_dataUS<- list_dataUS  %>% 
  mutate(treatfac=as.factor(treated),
         IDvar=as.factor(IDvar),
         homonat=as.factor(homonat),
         round.fac=as.factor(round))
listUS_A <- list_dataUS %>% 
  filter(round==1)
listUS_B <- list_dataUS %>% 
  filter(round==2)



###TABLE A.11###
##Pooled lists - Unadjusted##
pooledmodels <- list(
  "NL" = lm(value ~ treatfac, weights=weight, data = list_dataNL),
  
  "DE" = lm(value ~ treatfac, weights=weight, data = list_dataDE),
  
  "UK" = lm(value ~ treatfac, data = list_dataUK),
  
  "US" = lm(value ~ treatfac, data = list_dataUS))
modelsummary(pooledmodels,cluster = "IDvar", stars=TRUE,output = "tables/TableA11.tex")


###TABLE A.12###
##Individual lists - Unadjusted##

listmodels_noadj <- list(
  "NLA" = lm(value ~ treatfac,
             weights=weight, data = listNL_A),
  "NLB" = lm(value ~ treatfac, 
             weights=weight, data = listNL_B),
  "DEA" = lm(value ~ treatfac, 
             weights=weight, data = listDE_A),
  "DEB" = lm(value ~ treatfac, 
             weights=weight, data = listDE_B),
  "UKA" = lm(value ~ treatfac, data = listUK_A),
  "UKB" = lm(value ~ treatfac, data = listUK_B),
  "USA" = lm(value ~ treatfac, data = listUS_A),
  "USB" = lm(value ~ treatfac, data = listUS_B))
modelsummary(listmodels_noadj, stars=TRUE, output = "tables/TableA12.tex")


###TABLE A.13###
##Pooled lists - Adjusted##
pooledmodels <- list(
  "NL" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, 
            weights=weight, data = list_dataNL),
  
  "DE" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, 
            weights=weight, data = list_dataDE),
  
  "UK" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, data = list_dataUK),
  
  "US" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, data = list_dataUS))
modelsummary(pooledmodels, cluster ="IDvar", stars=TRUE, output = "tables/TableA13.tex")

###TABLE A.14###
##Individual lists - Adjusted##
listmodels <- list(
  "NLA" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
               respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, 
             weights=weight, data = listNL_A),
  "NLB" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
               respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, 
             weights=weight, data = listNL_B),
  
  "DEA" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
               respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, 
             weights=weight, data = listDE_A),
  "DEB" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
               respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, 
             weights=weight, data = listDE_B),
  
  "UKA" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
               respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, data = listUK_A),
  
  "UKB" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
               respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, data = listUK_B),
  
  "USA" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
               respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, data = listUS_A),
  
  "USB" = lm(value ~ treatfac + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
               respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality, data = listUS_B))
modelsummary(listmodels, stars=TRUE, output = "tables/TableA14.tex")


###TABLEA.15###
LGBTmodels <- list(
  "Netherlands"<- lm(value ~ treatfac*respondent_sexuality + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
                       respondent_ideology + affective_muslims + affective_lgb, 
                     weights=weight, data = list_dataNL),
  "Germany"<- lm(value ~ treatfac*respondent_sexuality + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
                   respondent_ideology + affective_muslims + affective_lgb, 
                 weights=weight, data = list_dataDE),
  "UK"<- lm(value ~ treatfac*respondent_sexuality + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb, 
            data = list_dataUK),
  "US"<- lm(value ~ treatfac*respondent_sexuality + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb, 
            data = list_dataUS))
modelsummary(LGBTmodels, cluster="IDvar", stars=TRUE, output="tables/TableA15.tex")



###FIGURE A3 & A4###
list_dataNL<- list_dataNL  %>% 
  mutate(respondent_sexuality=as.factor(respondent_sexuality))
list_dataDE<- list_dataDE  %>% 
  mutate(respondent_sexuality=as.factor(respondent_sexuality))
list_dataUK<- list_dataUK  %>% 
  mutate(respondent_sexuality=as.factor(respondent_sexuality))
list_dataUS<- list_dataUS  %>% 
  mutate(respondent_sexuality=as.factor(respondent_sexuality))

NL1X<- lm(value ~ treatfac*respondent_sexuality + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_ideology + affective_muslims + affective_lgb, 
          weights=weight, data = list_dataNL)
NLtidyX <- tidy(NL1X)
NLtidyX$cntry <- "Netherlands"
DE1X<- lm(value ~ treatfac*respondent_sexuality + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_ideology + affective_muslims + affective_lgb, 
          weights=weight, data = list_dataDE)
DEtidyX <- tidy(DE1X)
DEtidyX$cntry <- "Germany"
UK1X<- lm(value ~ treatfac*respondent_sexuality + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_ideology + affective_muslims + affective_lgb, 
          data = list_dataUK)
UKtidyX <- tidy(UK1X)
UKtidyX$cntry <- "UK"
US1X<- lm(value ~ treatfac*respondent_sexuality + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_ideology + affective_muslims + affective_lgb, 
          data = list_dataUS)
UStidyX <- tidy(US1X)
UStidyX$cntry <- "USA"

I2<- cat_plot(model = NL1X, modx = treatfac,
              pred = respondent_sexuality, robust=TRUE,
              interval.geom = "linerange",
              modx.labels = c("Control", "Treatment"))+
  scale_x_discrete(labels=c( "0"= "Cis-Hetero", "1" = "LGBTQ+"))+
  labs(title="Netherlands", y="Item-Count",
       x="Sexuality")+
  ylim(1.5, 3)+
  scale_colour_manual(values = c("hotpink", "pink"))+
  theme(legend.position = "none")

I1<- cat_plot(model = DE1X, modx = treatfac,
              pred = respondent_sexuality, robust=TRUE,
              interval.geom = "linerange",
              modx.labels = c("Control", "Treatment"))+
  scale_x_discrete(labels=c( "0"= "Cis-Hetero", "1" = "LGBTQ+"))+
  labs(title="Germany", y="Item-Count",
       x="Sexuality")+
  ylim(1.5, 3)+
  annotate(
    geom="text", x = 1, y = 2.1, size = 4, color = "hotpink", fontface=2,
    hjust=0,
    label = "Control")+
  annotate(
    geom="text", x = 1, y = 2.2, size = 4, color = "pink", fontface=2,
    hjust=0, label = "Treatment")+
  scale_colour_manual(values = c("hotpink", "pink"))+
  theme(legend.position = "none")

I3<- cat_plot(model = UK1X, modx = treatfac,
              pred = respondent_sexuality, robust=TRUE,
              interval.geom = "linerange",
              modx.labels = c("Control", "Treatment"))+
  scale_x_discrete(labels=c( "0"= "Cis-Hetero", "1" = "LGBTQ+"))+
  labs(title="UK", y="Item-Count",
       x="Sexuality")+
  ylim(1.5, 3)+
  scale_colour_manual(values = c("hotpink", "pink"))+
  theme(legend.position = "none")

I4<- cat_plot(model = US1X, modx = treatfac,
              pred = respondent_sexuality, robust=TRUE,
              interval.geom = "linerange",
              modx.labels = c("Control", "Treatment"))+
  scale_x_discrete(labels=c( "0"= "Cis-Hetero", "1" = "LGBTQ+"))+
  labs(title="US", y="Item-Count",
       x="Sexuality")+
  ylim(1.5, 3)+
  scale_colour_manual(values = c("hotpink", "pink"))+
  theme(legend.position = "none")

(I1+I2)/(I3+I4)
ggsave("FigureA3.png",
       path = here("figures"),
       dpi=600)


interact<- rbind(NLtidyX, DEtidyX, UKtidyX, UStidyX)

interact<-subset(interact, term %in% c("treatfac1:respondent_sexuality1")) 

interact$lower90 <- interact$estimate - (1.645 * interact$std.error)
interact$higher90 <- interact$estimate + (1.645 * interact$std.error)
interact$lower95 <- interact$estimate - (1.96 * interact$std.error)
interact$higher95 <- interact$estimate + (1.96 * interact$std.error)
interact$lower99 <- interact$estimate - (2.58 * interact$std.error)
interact$higher99 <- interact$estimate + (2.58 * interact$std.error)
interact$lower999 <- interact$estimate - (3.29 * interact$std.error)
interact$higher999 <- interact$estimate + (3.29 * interact$std.error)

ggplot(interact, aes(y = reorder(cntry, -estimate), x=estimate, color=cntry)) + 
  geom_vline(xintercept = 0, linetype="longdash", colour="black")+
  scale_color_manual(values = colors, guide = "none")+
  geom_segment(aes(x=lower90, xend=higher90, y=cntry, yend=cntry), size=3.2)+ 
  geom_segment(aes(x=lower95, xend=higher95, y=cntry, yend=cntry), size=2.5, alpha=.8)+ 
  geom_segment(aes(x=lower99, xend=higher99, y=cntry, yend=cntry), size=2, alpha=.7)+
  geom_segment(aes(x=lower999, xend=higher999, y=cntry, yend=cntry), size=1.5, alpha=.6)+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  xlim(-1,1)+
  geom_label(aes(label = cntry), nudge_y=.2)+
  geom_label(aes(label=sprintf("Diff. in ATE = %.2f", round(estimate, digits = 2))), nudge_y = -.15, size=2)+
  theme_bw()+
  labs(x="<<ATE lower than cis-heterosexuals -- ATE higher than cis-heterosexuals>>",
       title="(No) Treatment Heteogeneity Between Cis-Hetosexual & LGBTQ+ Respondents",
       caption= "Confidence intervals at 90%, 95%, 99%, and 99.9%")+
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        plot.title=element_text(color="black", face="bold"),
        axis.text.x=element_text(color="black"),
        axis.title.x=element_text(color="black", face="bold"),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())


ggsave("FigureA4.png",
       path = here("figures"),
       dpi=600)


###FIGUREA5###
NL1X<- lm(value ~ treatfac*respondent_ideology + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + affective_muslims + affective_lgb, 
          weights=weight, data = list_dataNL)


margNL<-
  NL1X %>%
  margins(at = list(respondent_ideology = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margNL$lower90 <- margNL$AME - (1.645 * margNL$SE)
margNL$higher90 <- margNL$AME + (1.645 * margNL$SE)
margNL$lower95 <- margNL$AME - (1.96 * margNL$SE)
margNL$higher95 <- margNL$AME + (1.96 * margNL$SE)
margNL$lower99 <- margNL$AME - (2.58 * margNL$SE)
margNL$higher99 <- margNL$AME + (2.58 * margNL$SE)
margNL$lower999 <- margNL$AME - (3.29 * margNL$SE)
margNL$higher999 <- margNL$AME + (3.29 * margNL$SE)


MP2 <-ggplot(margNL, aes(respondent_ideology, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataNL, aes(y = ..density..*1.5), fill = "#E16462FF", color = "#E16462FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=respondent_ideology, xend=respondent_ideology), size=3.2, color="#E16462FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=respondent_ideology, xend=respondent_ideology), size=2.5, alpha=.8, color="#E16462FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=respondent_ideology, xend=respondent_ideology), size=2, alpha=.7, color="#E16462FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=respondent_ideology, xend=respondent_ideology), size=1.5, alpha=.6, color="#E16462FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margNL, aes(label = round(AME, 2)), nudge_y = -0.35, color = "#E16462FF", size=2.5) +
  labs(title="Netherlands",
       x="Left-right placement",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank())


DE1X<- lm(value ~ treatfac*respondent_ideology + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + affective_muslims + affective_lgb, 
          weights=weight, data = list_dataDE)

margDE<-
  DE1X %>%
  margins(at = list(respondent_ideology = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margDE$lower90 <- margDE$AME - (1.645 * margDE$SE)
margDE$higher90 <- margDE$AME + (1.645 * margDE$SE)
margDE$lower95 <- margDE$AME - (1.96 * margDE$SE)
margDE$higher95 <- margDE$AME + (1.96 * margDE$SE)
margDE$lower99 <- margDE$AME - (2.58 * margDE$SE)
margDE$higher99 <- margDE$AME + (2.58 * margDE$SE)
margDE$lower999 <- margDE$AME - (3.29 * margDE$SE)
margDE$higher999 <- margDE$AME + (3.29 * margDE$SE)

MP1 <-ggplot(margDE, aes(respondent_ideology, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataDE, aes(y = ..density..*1.5), fill = "#FCA636FF", color = "#FCA636FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=respondent_ideology, xend=respondent_ideology), size=3.2, color="#FCA636FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=respondent_ideology, xend=respondent_ideology), size=2.5, alpha=.8, color="#FCA636FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=respondent_ideology, xend=respondent_ideology), size=2, alpha=.7, color="#FCA636FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=respondent_ideology, xend=respondent_ideology), size=1.5, alpha=.6, color="#FCA636FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margDE, aes(label = round(AME, 2)), nudge_y = -0.35, color = "#FCA636FF", size=2.5) +
  labs(title="Germany",
       x="Left-right placement",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x = element_blank(),
        axis.title.x = element_blank())


UK1X<- lm(value ~ treatfac*respondent_ideology + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + affective_muslims + affective_lgb,
          data = list_dataUK)

margUK<-
  UK1X %>%
  margins(at = list(respondent_ideology = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margUK$lower90 <- margUK$AME - (1.645 * margUK$SE)
margUK$higher90 <- margUK$AME + (1.645 * margUK$SE)
margUK$lower95 <- margUK$AME - (1.96 * margUK$SE)
margUK$higher95 <- margUK$AME + (1.96 * margUK$SE)
margUK$lower99 <- margUK$AME - (2.58 * margUK$SE)
margUK$higher99 <- margUK$AME + (2.58 * margUK$SE)
margUK$lower999 <- margUK$AME - (3.29 * margUK$SE)
margUK$higher999 <- margUK$AME + (3.29 * margUK$SE)


MP3 <-ggplot(margUK, aes(respondent_ideology, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataUK, aes(y = ..density..*1.5), fill = "#B12A90FF", color = "#B12A90FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=respondent_ideology, xend=respondent_ideology), size=3.2, color="#B12A90FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=respondent_ideology, xend=respondent_ideology), size=2.5, alpha=.8, color="#B12A90FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=respondent_ideology, xend=respondent_ideology), size=2, alpha=.7, color="#B12A90FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=respondent_ideology, xend=respondent_ideology), size=1.5, alpha=.6, color="#B12A90FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margUK, aes(label = round(AME, 2)), nudge_y = -0.3, color = "#B12A90FF", size=2.5) +
  labs(title="UK",
       x="Left-right placement",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5))


US1X<- lm(value ~ treatfac*respondent_ideology + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + affective_muslims + affective_lgb,
          data = list_dataUS)

margUS<-
  US1X %>%
  margins(at = list(respondent_ideology = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margUS$lower90 <- margUS$AME - (1.645 * margUS$SE)
margUS$higher90 <- margUS$AME + (1.645 * margUS$SE)
margUS$lower95 <- margUS$AME - (1.96 * margUS$SE)
margUS$higher95 <- margUS$AME + (1.96 * margUS$SE)
margUS$lower99 <- margUS$AME - (2.58 * margUS$SE)
margUS$higher99 <- margUS$AME + (2.58 * margUS$SE)
margUS$lower999 <- margUS$AME - (3.29 * margUS$SE)
margUS$higher999 <- margUS$AME + (3.29 * margUS$SE)


MP4 <-ggplot(margUS, aes(respondent_ideology, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataUS, aes(y = ..density..*1.5), fill = "#6A00A8FF", color = "#6A00A8FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=respondent_ideology, xend=respondent_ideology), size=3.2, color="#6A00A8FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=respondent_ideology, xend=respondent_ideology), size=2.5, alpha=.8, color="#6A00A8FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=respondent_ideology, xend=respondent_ideology), size=2, alpha=.7, color="#6A00A8FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=respondent_ideology, xend=respondent_ideology), size=1.5, alpha=.6, color="#6A00A8FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margUS, aes(label = round(AME, 2)), nudge_y = -0.3, color = "#6A00A8FF", size=2.5) +
  labs(title="USA",
       x="Left-right placement",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.y = element_blank(),
        axis.title.y = element_blank())
(MP1+MP2)/(MP3+MP4)+
  plot_annotation(caption= "Confidence intervals at 90%, 95%, 99%, and 99.9%")
ggsave("FigureA5.png", 
       path = here("figures"),
       dpi=1000)


###FIGURE A6###

NL1X<- lm(value ~ treatfac*affective_muslims + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + respondent_ideology + affective_lgb, 
          weights=weight, data = list_dataNL)

margNL<-
  NL1X %>%
  margins(at = list(affective_muslims = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margNL$lower90 <- margNL$AME - (1.645 * margNL$SE)
margNL$higher90 <- margNL$AME + (1.645 * margNL$SE)
margNL$lower95 <- margNL$AME - (1.96 * margNL$SE)
margNL$higher95 <- margNL$AME + (1.96 * margNL$SE)
margNL$lower99 <- margNL$AME - (2.58 * margNL$SE)
margNL$higher99 <- margNL$AME + (2.58 * margNL$SE)
margNL$lower999 <- margNL$AME - (3.29 * margNL$SE)
margNL$higher999 <- margNL$AME + (3.29 * margNL$SE)


MP2 <-ggplot(margNL, aes(affective_muslims, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataNL, aes(y = ..density..), fill = "#E16462FF", color = "#E16462FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=affective_muslims, xend=affective_muslims), size=3.2, color="#E16462FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=affective_muslims, xend=affective_muslims), size=2.5, alpha=.8, color="#E16462FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=affective_muslims, xend=affective_muslims), size=2, alpha=.7, color="#E16462FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=affective_muslims, xend=affective_muslims), size=1.5, alpha=.6, color="#E16462FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margNL, aes(label = round(AME, 2)), nudge_y = -0.35, color = "#E16462FF", size=2.5) +
  labs(title="Netherlands",
       x="Affect towards Muslims",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank())


DE1X<- lm(value ~ treatfac*affective_muslims + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + respondent_ideology + affective_lgb, 
          weights=weight, data = list_dataDE)

margDE<-
  DE1X %>%
  margins(at = list(affective_muslims = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margDE$lower90 <- margDE$AME - (1.645 * margDE$SE)
margDE$higher90 <- margDE$AME + (1.645 * margDE$SE)
margDE$lower95 <- margDE$AME - (1.96 * margDE$SE)
margDE$higher95 <- margDE$AME + (1.96 * margDE$SE)
margDE$lower99 <- margDE$AME - (2.58 * margDE$SE)
margDE$higher99 <- margDE$AME + (2.58 * margDE$SE)
margDE$lower999 <- margDE$AME - (3.29 * margDE$SE)
margDE$higher999 <- margDE$AME + (3.29 * margDE$SE)


MP1 <-ggplot(margDE, aes(affective_muslims, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataDE, aes(y = ..density..), fill = "#FCA636FF", color = "#FCA636FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=affective_muslims, xend=affective_muslims), size=3.2, color="#FCA636FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=affective_muslims, xend=affective_muslims), size=2.5, alpha=.8, color="#FCA636FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=affective_muslims, xend=affective_muslims), size=2, alpha=.7, color="#FCA636FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=affective_muslims, xend=affective_muslims), size=1.5, alpha=.6, color="#FCA636FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margDE, aes(label = round(AME, 2)), nudge_y = -0.35, color = "#FCA636FF", size=2.5) +
  labs(title="Germany",
       x="Affect towards Muslims",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x = element_blank(),
        axis.title.x = element_blank())


UK1X<- lm(value ~ treatfac*affective_muslims + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + respondent_ideology + affective_lgb,
          data = list_dataUK)

margUK<-
  UK1X %>%
  margins(at = list(affective_muslims = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margUK$lower90 <- margUK$AME - (1.645 * margUK$SE)
margUK$higher90 <- margUK$AME + (1.645 * margUK$SE)
margUK$lower95 <- margUK$AME - (1.96 * margUK$SE)
margUK$higher95 <- margUK$AME + (1.96 * margUK$SE)
margUK$lower99 <- margUK$AME - (2.58 * margUK$SE)
margUK$higher99 <- margUK$AME + (2.58 * margUK$SE)
margUK$lower999 <- margUK$AME - (3.29 * margUK$SE)
margUK$higher999 <- margUK$AME + (3.29 * margUK$SE)


MP3 <-ggplot(margUK, aes(affective_muslims, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataUK, aes(y = ..density..), fill = "#B12A90FF", color = "#B12A90FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=affective_muslims, xend=affective_muslims), size=3.2, color="#B12A90FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=affective_muslims, xend=affective_muslims), size=2.5, alpha=.8, color="#B12A90FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=affective_muslims, xend=affective_muslims), size=2, alpha=.7, color="#B12A90FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=affective_muslims, xend=affective_muslims), size=1.5, alpha=.6, color="#B12A90FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margUK, aes(label = round(AME, 2)), nudge_y = -0.3, color = "#B12A90FF", size=2.5) +
  labs(title="UK",
       x="Affect towards Muslims",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5))


US1X<- lm(value ~ treatfac*affective_muslims + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + respondent_ideology + affective_lgb,
          data = list_dataUS)

margUS<-
  US1X %>%
  margins(at = list(affective_muslims = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margUS$lower90 <- margUS$AME - (1.645 * margUS$SE)
margUS$higher90 <- margUS$AME + (1.645 * margUS$SE)
margUS$lower95 <- margUS$AME - (1.96 * margUS$SE)
margUS$higher95 <- margUS$AME + (1.96 * margUS$SE)
margUS$lower99 <- margUS$AME - (2.58 * margUS$SE)
margUS$higher99 <- margUS$AME + (2.58 * margUS$SE)
margUS$lower999 <- margUS$AME - (3.29 * margUS$SE)
margUS$higher999 <- margUS$AME + (3.29 * margUS$SE)


MP4 <-ggplot(margUS, aes(affective_muslims, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataUS, aes(y = ..density..), fill = "#6A00A8FF", color = "#6A00A8FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=affective_muslims, xend=affective_muslims), size=3.2, color="#6A00A8FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=affective_muslims, xend=affective_muslims), size=2.5, alpha=.8, color="#6A00A8FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=affective_muslims, xend=affective_muslims), size=2, alpha=.7, color="#6A00A8FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=affective_muslims, xend=affective_muslims), size=1.5, alpha=.6, color="#6A00A8FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margUS, aes(label = round(AME, 2)), nudge_y = -0.3, color = "#6A00A8FF", size=2.5) +
  labs(title="USA",
       x="Affect towards Muslims",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.y = element_blank(),
        axis.title.y = element_blank())

(MP1+MP2)/(MP3+MP4)+
  plot_annotation(caption= "Confidence intervals at 90%, 95%, 99%, and 99.9%")
ggsave("FigureA6.png", 
       path = here("figures"),
       dpi=1000)




###FIGURE A7###

NL1X<- lm(value ~ treatfac*affective_lgb + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + respondent_ideology + affective_muslims, 
          weights=weight, data = list_dataNL)

margNL<-
  NL1X %>%
  margins(at = list(affective_lgb = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margNL$lower90 <- margNL$AME - (1.645 * margNL$SE)
margNL$higher90 <- margNL$AME + (1.645 * margNL$SE)
margNL$lower95 <- margNL$AME - (1.96 * margNL$SE)
margNL$higher95 <- margNL$AME + (1.96 * margNL$SE)
margNL$lower99 <- margNL$AME - (2.58 * margNL$SE)
margNL$higher99 <- margNL$AME + (2.58 * margNL$SE)
margNL$lower999 <- margNL$AME - (3.29 * margNL$SE)
margNL$higher999 <- margNL$AME + (3.29 * margNL$SE)


MP2 <-ggplot(margNL, aes(affective_lgb, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataNL, aes(y = ..density..), fill = "#E16462FF", color = "#E16462FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=affective_lgb, xend=affective_lgb), size=3.2, color="#E16462FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=affective_lgb, xend=affective_lgb), size=2.5, alpha=.8, color="#E16462FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=affective_lgb, xend=affective_lgb), size=2, alpha=.7, color="#E16462FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=affective_lgb, xend=affective_lgb), size=1.5, alpha=.6, color="#E16462FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margNL, aes(label = round(AME, 2)), nudge_y = -0.35, color = "#E16462FF", size=2.5) +
  labs(title="Netherlands",
       x="Affect towards LGBT+",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank())


DE1X<- lm(value ~ treatfac*affective_lgb + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + respondent_ideology + affective_muslims,
          weights=weight, data = list_dataDE)

margDE<-
  DE1X %>%
  margins(at = list(affective_lgb = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margDE$lower90 <- margDE$AME - (1.645 * margDE$SE)
margDE$higher90 <- margDE$AME + (1.645 * margDE$SE)
margDE$lower95 <- margDE$AME - (1.96 * margDE$SE)
margDE$higher95 <- margDE$AME + (1.96 * margDE$SE)
margDE$lower99 <- margDE$AME - (2.58 * margDE$SE)
margDE$higher99 <- margDE$AME + (2.58 * margDE$SE)
margDE$lower999 <- margDE$AME - (3.29 * margDE$SE)
margDE$higher999 <- margDE$AME + (3.29 * margDE$SE)


MP1 <-ggplot(margDE, aes(affective_lgb, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataDE, aes(y = ..density..), fill = "#FCA636FF", color = "#FCA636FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=affective_lgb, xend=affective_lgb), size=3.2, color="#FCA636FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=affective_lgb, xend=affective_lgb), size=2.5, alpha=.8, color="#FCA636FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=affective_lgb, xend=affective_lgb), size=2, alpha=.7, color="#FCA636FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=affective_lgb, xend=affective_lgb), size=1.5, alpha=.6, color="#FCA636FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margDE, aes(label = round(AME, 2)), nudge_y = -0.35, color = "#FCA636FF", size=2.5) +
  labs(title="Germany",
       x="Affect towards LGBT+",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x = element_blank(),
        axis.title.x = element_blank())


UK1X<- lm(value ~ treatfac*affective_lgb + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + respondent_ideology + affective_muslims,
          data = list_dataUK)

margUK<-
  UK1X %>%
  margins(at = list(affective_lgb = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margUK$lower90 <- margUK$AME - (1.645 * margUK$SE)
margUK$higher90 <- margUK$AME + (1.645 * margUK$SE)
margUK$lower95 <- margUK$AME - (1.96 * margUK$SE)
margUK$higher95 <- margUK$AME + (1.96 * margUK$SE)
margUK$lower99 <- margUK$AME - (2.58 * margUK$SE)
margUK$higher99 <- margUK$AME + (2.58 * margUK$SE)
margUK$lower999 <- margUK$AME - (3.29 * margUK$SE)
margUK$higher999 <- margUK$AME + (3.29 * margUK$SE)


MP3 <-ggplot(margUK, aes(affective_lgb, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataUK, aes(y = ..density..), fill = "#B12A90FF", color = "#B12A90FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=affective_lgb, xend=affective_lgb), size=3.2, color="#B12A90FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=affective_lgb, xend=affective_lgb), size=2.5, alpha=.8, color="#B12A90FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=affective_lgb, xend=affective_lgb), size=2, alpha=.7, color="#B12A90FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=affective_lgb, xend=affective_lgb), size=1.5, alpha=.6, color="#B12A90FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margUK, aes(label = round(AME, 2)), nudge_y = -0.3, color = "#B12A90FF", size=2.5) +
  labs(title="UK",
       x="Affect towards LGBT+",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5))


US1X<- lm(value ~ treatfac*affective_lgb + respondent_age + respondent_age2 + respondent_gender + respondent_degree +
            respondent_sexuality + respondent_ideology + affective_muslims,
          data = list_dataUS)

margUS<-
  US1X %>%
  margins(at = list(affective_lgb = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treatfac1")

margUS$lower90 <- margUS$AME - (1.645 * margUS$SE)
margUS$higher90 <- margUS$AME + (1.645 * margUS$SE)
margUS$lower95 <- margUS$AME - (1.96 * margUS$SE)
margUS$higher95 <- margUS$AME + (1.96 * margUS$SE)
margUS$lower99 <- margUS$AME - (2.58 * margUS$SE)
margUS$higher99 <- margUS$AME + (2.58 * margUS$SE)
margUS$lower999 <- margUS$AME - (3.29 * margUS$SE)
margUS$higher999 <- margUS$AME + (3.29 * margUS$SE)


MP4 <-ggplot(margUS, aes(affective_lgb, AME)) +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.1, 1)) +
  geom_density(data = list_dataUS, aes(y = ..density..), fill = "#6A00A8FF", color = "#6A00A8FF", alpha = 0.5) +  # Add the density plot
  geom_segment(aes(y=lower90, yend=higher90, x=affective_lgb, xend=affective_lgb), size=3.2, color="#6A00A8FF")+ 
  geom_segment(aes(y=lower95, yend=higher95, x=affective_lgb, xend=affective_lgb), size=2.5, alpha=.8, color="#6A00A8FF")+ 
  geom_segment(aes(y=lower99, yend=higher99, x=affective_lgb, xend=affective_lgb), size=2, alpha=.7, color="#6A00A8FF")+
  geom_segment(aes(y=lower999, yend=higher999, x=affective_lgb, xend=affective_lgb), size=1.5, alpha=.6, color="#6A00A8FF")+
  geom_point(fill="white", color="grey67", size=2.5,  pch=21)+
  geom_hline(yintercept = 0, linetype="longdash", colour="black")+
  geom_label(data = margUS, aes(label = round(AME, 2)), nudge_y = -0.3, color = "#6A00A8FF", size=2.5) +
  labs(title="USA",
       x="Affect towards LGBT+",
       y="Conditional ATE") +
  theme_minimal()+
  theme(plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.y = element_blank(),
        axis.title.y = element_blank())

(MP1+MP2)/(MP3+MP4)+
  plot_annotation(caption= "Confidence intervals at 90%, 95%, 99%, and 99.9%")
ggsave("FigureA7.png",
       path = here("figures"),
       dpi=1000)


###TABLE A16###
dt_DE<- ict.test(list_dataDE$value, list_dataDE$treated, J = 4, alpha = 0.05,
                 n.draws = 250000, gms = TRUE,  pi.table = TRUE)
coeffsDE<- dt_DE$pi.table
latex_tabDE <- data.frame(Coefficient = coeffsDE)
xtable(latex_tabDE)
xtable::xtable(latex_tabDE,
               digits = 3) -> esttab_latex_tabDE

xtable::print.xtable(esttab_latex_tabDE, 
                     include.rownames = T,
                     type = "latex",
                     file = "tables/TableA16.tex")


###TABLE A17###
dt_NL<- ict.test(list_dataNL$value, list_dataNL$treated, J = 4, alpha = 0.05,
                 n.draws = 250000, gms = TRUE,  pi.table = TRUE)
coeffsNL<- dt_NL$pi.table
latex_tabNL <- data.frame(Coefficient = coeffsNL)
xtable(latex_tabNL)
xtable::xtable(latex_tabNL,
               digits = 3) -> esttab_latex_tabNL

xtable::print.xtable(esttab_latex_tabNL, 
                     include.rownames = T,
                     type = "latex",
                     file = "tables/TableA17.tex")

###TABLE A18###
dt_UK<- ict.test(list_dataUK$value, list_dataUK$treated, J = 4, alpha = 0.05,
                 n.draws = 250000, gms = TRUE,  pi.table = TRUE)
coeffsUK<- dt_UK$pi.table
latex_tabUK <- data.frame(Coefficient = coeffsUK)
xtable(latex_tabUK)
xtable::xtable(latex_tabUK,
               digits = 3) -> esttab_latex_tabUK

xtable::print.xtable(esttab_latex_tabUK, 
                     include.rownames = T,
                     type = "latex",
                     file = "tables/TableA18.tex")

###TABLE A19###
dt_US<- ict.test(list_dataUS$value, list_dataUS$treated, J = 4, alpha = 0.05,
                 n.draws = 250000, gms = TRUE,  pi.table = TRUE)
coeffsUS<- dt_US$pi.table
latex_tabUS <- data.frame(Coefficient = coeffsUS)
xtable(latex_tabUS)
xtable::xtable(latex_tabUS,
               digits = 3) -> esttab_latex_tabUS

xtable::print.xtable(esttab_latex_tabUS, 
                     include.rownames = T,
                     type = "latex",
                     file = "tables/TableA19.tex")
###TABLE A20###
DD <- list(
  "NL" = lm_robust(value ~ treatfac*round.fac, clusters = IDvar, weights=weight,
                   data = list_dataNL),
  "DE" = lm_robust(value ~ treatfac*round.fac, clusters = IDvar, weights=weight,
                   data = list_dataDE),
  "UK" = lm_robust(value ~ treatfac*round.fac, clusters = IDvar,
                   data = list_dataUK),
  "US" = lm_robust(value ~ treatfac*round.fac, clusters = IDvar,
                   data = list_dataUS))

modelsummary(DD, star=TRUE, cluster ="IDvar",
             output='tables/TableA20.tex')


###FIGURE A12###

list_dataNL$cntry <- "Netherlands"
list_dataDE$cntry <- "Germany"
list_dataUK$cntry <- "UK"
list_dataUS$cntry <- "US"

list_df <- list(list_dataNL, list_dataDE, list_dataUK, list_dataUS)
fulldf <- bind_rows(list_df)

perm_controls = c(
  'Gender' = 'respondent_gender',
  'Age' = 'respondent_age',
  'LGBT+' = 'respondent_sexuality',
  'Education' = 'respondent_degree',
  'Ideology' = 'respondent_ideology',
  'Affect towards Muslims' = 'affective_muslims',
  'Affect towards LGBT+' = 'affective_lgb'
)

perm_fe_controls = c(
  'Country FE' = 'cntry',
  'List FE' = 'round'
)
perm_fe_controls2 = c(
  'List FE' = 'round'
)

p1<- stability_plot(data = fulldf, 
                    lhs = 'value', 
                    rhs = 'treatfac', 
                    perm = perm_controls,
                    perm_fe = perm_fe_controls,
                    sort = 'desc-by-fe',
                    error_geom = 'ribbon', 
                    coef_ylim = c(.45, .75),  # change the endpoints of the y-axis
                    rel_height = 0.6,
                    run_to = 6
)


replacement_coef_panel <- p1[[1]] + 
  labs(title="Combined sample") +
  scale_color_manual(values = c("#673B90"))+
  theme(plot.title = element_text(face="bold"))

replacement_control_panel = p1[[2]] + 
  scale_fill_manual(values = c("#FFFFFF", "#673B90"))

p1X<- combine_plots(replacement_coef_panel, 
                    replacement_control_panel, 
                    rel_height = 0.6)
ggsave("FigureA12.png", 
       path = here("figures"),
       dpi=1000, height=12, width=18, unit="cm")



###FIGURE A13
p1<- stability_plot(data = list_dataDE, 
                    lhs = 'value', 
                    rhs = 'treatfac', 
                    perm = perm_controls,
                    perm_fe = perm_fe_controls2,
                    sort = 'desc-by-fe',
                    error_geom = 'ribbon', 
                    coef_ylim = c(0.45, .75),  # change the endpoints of the y-axis
                    rel_height = 0.6,
                    run_to = 6
)


replacement_coef_panel <- p1[[1]] + 
  labs(title="Germany") +
  scale_color_manual(values = c("#673B90"))+
  theme(plot.title = element_text(face="bold"))

replacement_control_panel = p1[[2]] + 
  scale_fill_manual(values = c("#FFFFFF", "#673B90"))

p1X<- combine_plots(replacement_coef_panel, 
                    replacement_control_panel, 
                    rel_height = 0.6)

ggsave("FigureA13_panel1.png", 
       path = here("figures"),dpi=1000, height=12, width=18, unit="cm")

p1<- stability_plot(data = list_dataNL, 
                    lhs = 'value', 
                    rhs = 'treatfac', 
                    perm = perm_controls,
                    perm_fe = perm_fe_controls2,
                    sort = 'desc-by-fe',
                    error_geom = 'ribbon', 
                    coef_ylim = c(0.45, .75),  # change the endpoints of the y-axis
                    rel_height = 0.6,
                    run_to = 6
)


replacement_coef_panel <- p1[[1]] + 
  labs(title="Netherlands") +
  scale_color_manual(values = c("#673B90"))+
  theme(plot.title = element_text(face="bold"))

replacement_control_panel = p1[[2]] + 
  scale_fill_manual(values = c("#FFFFFF", "#673B90"))

p1X<- combine_plots(replacement_coef_panel, 
                    replacement_control_panel, 
                    rel_height = 0.6)

ggsave("FigureA13_panel2.png", 
       path = here("figures"),
       dpi=1000, height=12, width=18, unit="cm")

p1<- stability_plot(data = list_dataUK, 
                    lhs = 'value', 
                    rhs = 'treatfac', 
                    perm = perm_controls,
                    perm_fe = perm_fe_controls2,
                    sort = 'desc-by-fe',
                    error_geom = 'ribbon', 
                    coef_ylim = c(0.45, .75),  # change the endpoints of the y-axis
                    rel_height = 0.6,
                    run_to = 6
)


replacement_coef_panel <- p1[[1]] + 
  labs(title="UK") +
  scale_color_manual(values = c("#673B90"))+
  theme(plot.title = element_text(face="bold"))

replacement_control_panel = p1[[2]] + 
  scale_fill_manual(values = c("#FFFFFF", "#673B90"))

p1X<- combine_plots(replacement_coef_panel, 
                    replacement_control_panel, 
                    rel_height = 0.6)

ggsave("FigureA13_panel3.png", 
       path = here("figures"),
       dpi=1000, height=12, width=18, unit="cm")

p1<- stability_plot(data = list_dataUS, 
                    lhs = 'value', 
                    rhs = 'treatfac', 
                    perm = perm_controls,
                    perm_fe = perm_fe_controls2,
                    sort = 'desc-by-fe',
                    error_geom = 'ribbon', 
                    coef_ylim = c(0.45, .75),  # change the endpoints of the y-axis
                    rel_height = 0.6,
                    run_to = 6
)


replacement_coef_panel <- p1[[1]] + 
  labs(title="US") +
  scale_color_manual(values = c("#673B90"))+
  theme(plot.title = element_text(face="bold"))

replacement_control_panel = p1[[2]] + 
  scale_fill_manual(values = c("#FFFFFF", "#673B90"))

p1X<- combine_plots(replacement_coef_panel, 
                    replacement_control_panel, 
                    rel_height = 0.6)

ggsave("FigureA13_panel4.png",
       path = here("figures"),
       dpi=1000, height=12, width=18, unit="cm")


###FIGURE A1###
fulldf<- fulldf  %>% 
  mutate(homonat=as.factor(homonat))
data_summary <- fulldf %>%
  group_by(cntry, homonat) %>%
  summarize(count = n()) %>%
  group_by(cntry) %>%
  mutate(percentage = count / sum(count))

ggplot(data_summary, aes(y = cntry, x = percentage, fill = homonat)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("0" = "hotpink", "1" = "gray60", "2" = "gray70", "3" = "gray50"),
                    labels = c("0" = "Sexually Modern Nativists",
                               "1" = "Sexually Conservative Nativists",
                               "2" = "Sexually Conservative Cosmopolitans",
                               "3" = "Sexually Modern Cosmopolitans"))+
  labs(x = "Percentage", title = "Proportion of Sexually Modern Nativists") +
  theme_minimal()+
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        axis.text.y = element_text(face="bold"),
        plot.title = element_text(face="bold"),
        axis.title.y = element_blank(),
        axis.title.x = element_blank())+
  guides(fill = guide_legend(nrow = 2))

ggsave("FigureA1.png", 
       path = here("figures"),
       dpi=600)

